rm(list=ls())

# Load packages
library(haven)
library(tidyverse)
library(data.table)
library(reshape2)
library(sf)
library(tmap)
library(sjlabelled)

source("code/config.R")
source(paste0(CODE,"parameters.R"))=

# Load map
mpsz <- st_read(dsn = paste0(DATA,"master-plan-2014-subzone-boundary-web-shp"), 
                layer = "MP14_SUBZONE_WEB_PL")

# Merge data
df <- read_dta(paste(DATA, "panelODByWeekAdultsOnly.dta", sep=""))
df <- df %>% filter(period=="before") %>% select(-c(period,dt,id))
df_main2 <- df %>% group_by(pHome) %>% summarise(share_residents= sum(totwork, na.rm=T),share_residents_w= sum(totworkw, na.rm=T)) %>% mutate(share_residents= share_residents/sum(share_residents, na.rm=T),share_residents_w= share_residents_w/sum(share_residents_w, na.rm=T))
df_main2$`Share High Income`<- df_main2$share_residents*total_pop / (df_main2$share_residents*total_pop + df_main2$share_residents_w*total_pop_w)
df_main2$name<-get_labels(df_main2$pHome)
df_main2<- df_main2 %>% select(name,`High Income Share2`)
mpsz_out <- left_join(mpsz_out, df_main2, 
                      by = c("PLN_AREA_N" = "name"))

DTLStationCoords <- read.csv("data/DTLStationCoords.csv", header = TRUE, stringsAsFactors = FALSE)
stations1 <- st_as_sf(DTLStationCoords, coords = c('LON', 'LAT'), crs=4326)

dest_data$PLN_AREA_N<-dest_data$sub_area_d
mpsz<-left_join(mpsz,dest_data)

# Create Figure 1a
qtm(mpsz_out, fill="Share High Income")+tm_shape(stations1)+tm_bubbles(col = "blue", scale = .3) 
